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,— I In this paper we present an extended version of Hilbert-Huang transform, namely arbitrary-order 

Hilbert spectral analysis, to characterize the scale-invariant properties of a time series directly in an 
amplitude-frequency space. We first show numerically that due to a nonlinear distortion, traditional 
^\ methods require high-order harmonic components to represent nonlinear processes, except for the 

i-H Hilbert-based method. This will lead to an artificial energy flux from the low-frequency (large scale) 

to the high-frequency (small scale) part. Thus the power law, if it exists, is contaminated. We then 
compare the Hilbert method with structure functions (SF), detrended fluctuation analysis (DFA), 
and wavelet leader (WL) by analyzing fractional Brownian motion and synthesized multifractal 
time series. For the former simulation, we find that all methods provide comparable results. For 
the latter simulation, we perform simulations with an intermittent parameter /j, — 0.15. We find 
that the SF underestimates scaling exponent when q > 3. The Hilbert method provides a slight 
underestimation when q > 5. However, both DFA and WL overestimate the scaling exponents when 
q > 5. It seems that Hilbert and DFA methods provide better singularity spectra than SF and WL. 
We finally apply all methods to a passive scalar (temperature) data obtained from a jet experiment 
with a Taylor' s microscale Reynolds number Reiambda ~ 250. Due to the presence of strong 
ramp-cliff structures, the SF fails to detect the power law behavior. For the traditional method, 
the ramp-cliff structure causes a serious artificial energy flux from the low-frequency (large scale) 
to the high-frequency (small scale) part. Thus DFA and WL underestimate the scaling exponents. 
However, the Hilbert method provides scaling exponents £e(q) quite close to the one for longitudinal 
velocity, indicating a less intermittent passive scalar field than what was believed before. 
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16 PACS numbers: 94.05.Lk, 05.45.Tp, 02.50.Fz 

VO 17 I. INTRODUCTION 32 increments and characterizes the scale invariance in the 

C*i 33 physical domain. 

is Multifractal properties have been found in many fields, 34 Other widely used methods are wavelet-based meth- 

19 such as turbulence [THS], rainfall [1H2], financial time se- 35 ods, e.g. wavelet transform modulus maxima (WTMM), 

20 ries [Mil], physiology [12], etc. Conventionally, multi- 36 wavelet leader (WL), or gradient modulus wavelet pro- 

21 fractal properties of such time series are characterized 37 jection (GMWP), to extract the scaling exponents from 
^ 22 by the scaling exponents ((q), which are extracted by 38 a scaling time series [TM5U] . However, as we will show in 

structure function (SF) analysis: AV* — ft( q \ where 39 this paper, the wavelets share the same drawback with 

AVe = \V(x + £)- V(x)\ are the increment for scale sep- m Fourier transform, which requires high-order harmonic 

'd 25 aration £, and ((q) is a nonlinear function [2d 02]. The 41 components to represent nonlinear processes. 

26 function C(q) is linear for monoscaling processes and non- e 1 , ., , , -, , 

.. ci... 1 ,tt 1 .42 borne of us nave proposed recently a new methodology, 

27 linear for multifractal processes. We may also mention , ... , „.,, , , , , . m „ u 

, , . . n . , . ,__. . . J rr-n rr-^ , 43 namely arbitrary-order Hubert spectral analysis (HbA), 

28 the detrended fluctuation analysis (DFA) IMHT6I or the , J -I, , • , .. ,. ., • 

1; ._ , , . . _ ,. v ' } — ii=r , . , 44 to characterize the scale invariant properties directly in 

29 multifractal detrended fluctuation analysis 1171 , which ,., , c ptthtpti tj. • j j 

. i t r t • — • i 45 amplitude-frequency space \6l oi . It is an extended 

30 are sometimes also employed lor scaling time series anal- . c ., TT . 1U , TT c /,,„„< ,. , 

mi t^t^a i f • • -i mt? .... 46 version of the Hilbert- Huang transform (HH1), which 

31 ysis. Ine Ur A method is similar to bl s since it involves ., • • , , i. v , , ■, <■ ,■ r jr\ • 

J 47 provides a joint probability density function (pdt) in an 

48 amplitude-frequency space [31ij 33]- We have applied part 

49 of this new methodology to several different time series 

* yongxianghuang@gmail.com 50 to show its efficiency and validity: turbulence experimen- 

t |trancois.schmitt<guniv-iiiierrr] 51 tal database |3T], synthesized fractional Brownian mo- 
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tion (ffim) time series [32] , surf zone marine turbulence 
data [34], and river flow discharge data [35]. In this pa- 
per we consider in length and precisely this new method 
and its validation and calibration. We first introduce 
this new methodology in detail. We then validate and 
calibrate it by analyzing a synthesized multifractal time 
series. We finally consider a passive scalar (temperature) 
data set with strong ramp-cliff structures. Due to the 
presence of ramp-cliff structures, the classical SF analy- 
sis fails to detect the power law behavior [36] . Addition- 
ally, for traditional methods, such as Fourier transform, 
wavelet transform, high-order harmonics is required to 
represent these structures and leads to an artificial en- 
ergy flux from the large-scale part (low frequency) to the 
small-scale part (high frequency) [31)] . 

This paper is organized as follows. We present the def- 
inition of the arbitrary-order Hilbert spectral analysis in 
Sec. [TT] The classical structure function analysis, multi- 
fractal detrended fluctuation analysis and wavelet leader 
are also presented in this section. We then consider a 
nonlinear effect by using the classical Duffing equation 
to show the artificial high-order harmonic components 
required by the classical methods, for example, Fourier 
transform, and wavelet transform, in Sec. |HI| In Sec. |IV| 
we perform a comparison study of the arbitrary-order 
HSA with other methods by analyzing a fBm simula- 
tion and a synthesized multifractal time series. We then 
present in Sec. [V]an analysis of real temperature data ob- 
tained from a jet experiment. We finally draw the main 
conclusions in Sec. [Vll 



82 II. METHODOLOGY 

83 A. Arbitrary-order Hilbert spectral analysis 

84 1. Empirical mode decomposition 

as The most innovative part of the Hilbert-Huang trans- 

86 form is the so-called empirical mode decomposition 

87 (EMD). In the real world most of the signals are multi- 

88 components, which means that different scales can co- 

89 exist simultaneously [37H35] . This may be considered 

90 as fast oscillations superposed to slower ones at a local 

91 level [301 SI]. Meanwhile, for decomposition methods, a 

92 characteristic scale (CS) is always defined implicitly or 

93 explicitly before the decomposition. For example, the 

94 CS of the classical Fourier analysis is a period of sine 

95 wave. The CS of wavelet transform is the shape of the 

96 mother wavelet [33 . In the present method, the CS is 

97 defined as the distance between two successive maxima 

98 (respectively minima) points. Then the so-called intrin- 

99 sic mode functions (IMF) are proposed to represent each 
ioo mono-component signal. An IMF satisfies the following 
ioi two conditions: (i) the difference between the number of 

102 local extrema and the number of zero-crossings must be 

103 zero or one; (ii) the running mean value of the envelope 

104 defined by the local maxima and the envelope defined by 



105 the local minima is zero [35] I3H] • 

106 A subpart of the EMD algorithm, called "sifting pro- 

107 cess," is then designed to decompose a given signal into 
loa several IMF modes [58H4T)] . The first step of the sifting 
109 process is to identify all the local maxima (respectively, 
no minima) points for a given time series x(t). Once all the 
in local extrema points are identified, the upper envelope 
112 e max (£) an d the lower envelope e m { n (t) are constructed, 
in respectively, for the local maxima and minima points by 
u4 using a cubic spline algorithm. The mean between these 
us two envelopes is defined as 



mi(t) 



c(*) + e min (i) 



Thus the first component is estimated by 
h^t) = x(t) -mi(t) 



(1) 



(2) 



n7 Ideally, ft-i(t) should be an IMF as expected. However, 
us h\{t) may not satisfy the above-mentioned conditions to 
n9 be an IMF. The function hi{t) is then taken as a new 

120 time series and this sifting process is repeated j times, 

121 until h\j{t) is an IMF 

hij{t) = h 1{j _ 1} (t) - mij(t) (3) 

122 The first IMF component C\ (t) is then written as 

C l (t)=h lj (t) (4) 

123 and the residual fi(t) as 



n(t) = x(t) - d(t) 



(5) 



from the data x(t). The sifting procedure is then re- 
peated on the residual until r n (t) becomes a monotonic 
function or at most has one local extreme point, which 
means that no more IMF can be extracted from r n {t). 
There are finally n — 1 IMF modes with one residual 
r n (t). The original signal x(t) is rewritten at the end of 
the process as 



X (t) = j2 c *® 



,(*) 



(6) 
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131 To guarantee that the IMF modes retain enough physical 

132 sense, a certain stopping criterion has to be introduced to 

133 stop the sifting process properly. Different types of stop- 

134 ping criteria have been introduced by several authors [38 - 

135 HOI EH S3] • The first stopping criterion is a Cauchy-type 

136 convergence criterion. We introduce the standard devi- 

137 ation (SD), defined for two successive sifting processes 

138 as 



SD 






-0 \hi(j- 



,(t)-/li(t)l s 
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139 If a calculated SD is smaller than a given value, then 
uo the sifting stops, and gives an IMF. A typical value is 



141 0.2 ~ 0.3, proposed based on Huang et al.'s experiences 

142 [38j [39] . Another widely used criterion is based on three 

143 thresholds a, 9\, and 02, which are designed to guarantee 

144 globally small fluctuations meanwhile taking into account 

145 locally large excursions |40j . The mode amplitude and 
«6 evaluation function are given as 



a(t) = 



^max v^ ) ^min \t) 



and 



a(t) = |m(*)/o(*)| 



(8a) 



(8b) 
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where P means the Cauchy principle value [371 I2H GEH 
1ST] . For each mode function Ci{t), one can then construct 
the analytical signal [371 |ST], Q(t), as 



where 



C i (t) = C i (t)+jC i (t) = A i (t)J e *W 



Mt) = [c i (t) 2 + cf(t)} 1 / 2 

9i (t) = arctan ( ^44 J 



(10) 



(11) 



Hence the instantaneous frequency is defined as 

1 M t (t) 



2?r di 



(12) 



The original signal is finally represented [excluding the 
residual r n (t)] as 



n 



N 



x{t) = RJ2 Mty eM = R^^iWe 



jf Ui (t)dt 



(13) 



«8 So that the sifting is iterated until a(t) < 8\ for some 

149 prescribed fraction 1 — a of the total duration, while 

150 a{t) < 82 for the remaining fraction. The typical values 

151 proposed by Rilling et al. [JD] are a ss 0.05, 0\ ~ 0.05 

152 and 02 ~ 10 01, respectively based on their experience. In 

153 our practice, if one of these criteria is satisfied, then the 

154 sifting process will stop. We also set a maximal iteration 

155 number (e.g., 300) to avoid over-decomposing the time 

156 series. 

157 The above-described EMD algorithm performs the de- 

158 composition on a very local level in the physical do- 

159 main without a priori basis. This means that the 
wo present decomposition is a posteriori: The basis is in- 
i6i duced by the data itself [35J I3S1 EI]. It is thus a 

162 scale-based decomposition. Since its introduction, this 

163 method has attracted large interests in various research 

164 fields: waves [3H 031 25], biological applications [3SJ- 

165 05] , financial studies [3H] , meteorology and climate stud- 

166 ies [35J [SDHS1], mechanical engineering [55] [5SJ, acous- 

167 tics [57] , aquatic environment [58] , and turbulence [31] , to 
ies quote a few. More detail about the EMD algorithm can 
169 be found in several methodological papers 



where "R" means real part. The above procedure is the 
classical Hilbert spectral analysis [37J |HI] • The combina- 
tion of EMD and HSA is thus called Hilbert-Huang trans- 
form by some authors |43] . The Hilbert-Huang transform 
can be taken as a generalization of the Fourier trans- 
form, since it allows frequency modulation and ampli- 
tude modulation simultaneously. The Hilbert spectrum, 
H(uj,t) = A 2 (uj,t), is designed to represent the energy 
in a time- frequency representation [60] . We further can 
define the Hilbert marginal spectrum as 



h(u) 



+OO 



H(uj,t)dt 



(14) 



170 2. Hilbert spectral analysis 

171 After having extracted the IMF modes, one can ap- 

172 ply the associated Hilbert spectral analysis to each com- 

173 ponent Ci in order to extract the energy time-frequency 

174 information from the data [351 GM ED] • The Hilbert trans- 

175 form of a function C(t) is written as 



This is similar with the Fourier spectrum, and can be 
interpreted as the energy associated with each frequency. 
However, we underline the fact that the definition of fre- 
quency here is different from the definition in the Fourier 
frame. Thus the interpretation of the Hilbert marginal 
spectrum should be given more caution J35J |3j5] • 



3. Arbitrary- order Hilbert spectral analysis 

We can also define the joint pdf p(u>, A) of the instanta- 
neous frequency uj and the amplitude A for each of these 
IMF modes [3TH331 I5D] . The Hilbert marginal spectrum 
Eq. ( 14 ) is then rewritten as 



h(u) 



+ OO 



p{iv,A)A 2 dA 



(15) 







204 The above definition is no more than the second-order 

205 statistical moment. This constatation has led some of 



206 



us to recently generalize this approach to arbitrary-order 
moment q > [3TH33] 



(16) 



/•+00 

C q (u)= / p{u 1 A)A q c\A 
Jo 

208 In case of scale invariance, we have 
£ g (w)~«-««> 



(17) 



209 in which £(g) is the Hilbert-based scaling exponent func- 

210 tion. Due to the integration operator, £ (q) — 1 can be 

211 associated with ((q) from SF analysis [31"ll33]. 

212 A limitation of the Hilbert-based method we proposed 

213 here is that it lacks the ability to consider q < [Q2 . In 



214 other words, similarly with the SF analysis, it has no res- 

215 olution on the right part of the singularity spectrum. The 

216 main drawback of the Hilbert-based method is its absence 

217 of solid theoretical ground, since the EMD part is almost 

218 empirical |43| . It has been found experimentally that 

219 the method, especially for the HSA, is statistically stable 

220 with different stopping criteria |42j . Recently, Flandrin 

221 et al. have obtained new theoretical results on the EMD 

222 method [4T1 fS"9l IB"3Tl6"5] . However, more theoretical work 

223 is still needed to fully mathematically understand this 

224 method. 



225 B. Structure function analysis 

226 The conventional way to extract scaling exponents is 

227 the classical SF analysis, which has been proposed in the 

228 field of turbulence and is now quite classical for intermit- 

229 tency studies [T3] . The qth order SF is written as 



S q (£) = (\Axt(t)\«) ~ £SM 



(18) 



Y(j)=J2(x(i)-x), j = l, 



■N 



(19) 



260 where x is the mean value of x. We then divide it into 
26i M n segments of length n (n < N) starting from both the 

262 beginning and the end of the time series. Each segment 

263 v has its own local trend that can be approximated by 

264 fitting a pth-order polynomial P% which is removed from 

265 the data. The variances for all the segments v and for all 

266 segment lengths n are then calculated by 



n 

F 2 (v,n) = -£{y[(„- i) n + j] -P*(j)} 2 



(20) 



i=i 



The gth-order fluctuation function is then defined as 



F q (n) = 



' , 2M„ 

2Mn lhi 



i/'-/ 



[F 2 (v,n)}' 



(21) 



For discussion convenience, we redefine the gth-order 
fluctuation function as 



230 where Axi(t) = x(t + I) — x(t) and £ is the time sepa- 

231 ration. The scaling exponent £(g) characterizes the fluc- 

232 tuation statistic at all scales; it is linear for monofractal 

233 processes such as fractional Brownian motion, and non- 
234 linear and concave (as a second Laplace characteristic 

235 function) for multifractal processes |66| . This approach 

236 has been widely used in turbulent research [U [2J |T3] and 

237 also other research fields [9l|67l|68]- However, the incre- 

238 ment operation acts a filter and thus SF characterizes the 

239 scale-invariant properties in an indirect way; see detailed 

240 discussion in Refs. (33J [3B] . 

241 As we have shown elsewhere, the increment operation 

242 in SF acts a filter and is a global operation. It thus 

243 measures the scale invariant property in an indirect way. 

244 It is also found that it is strongly influenced by energetic 

245 large scale structures [33, 36 . Therefore the SF analysis 

246 is not suitable for those data which possess energetic large 

247 scale structures. We will show an example of passive 

248 scalar turbulence data with strong ramp-cliff structures 

249 in Sec. [VJ More discussion can be found in Refs. [33| 136] . 

250 C. Multifractal detrended fluctuation analysis 

251 DFA was first introduced by Peng et al. [T3] to study 

252 the scaling properties of DNA sequence, in which only 

253 the second-order moment q = 2 was considered. Later 

254 this was generalized into a multifractal version by con- 

255 sidering the arbitrary order q, namely multifractal de- 

256 trended fluctuation analysis (MFDFA) [17l [69]. It then 

257 became a more common technique for scaling data anal- 

258 ysis [TlHrTl |6T)HT5] . For a given discrete time series x(i), 

259 i = 1 • • ■ A, we first estimate its cumulative function 



T q {n) = F"{nf 



(22) 



270 In case of scale invariance, we have power law scaling 

271 within a significant range of n 



T q {n) 



n 



h(q) 



(23) 



272 in which h{q) is the corresponding scaling exponent func- 

273 tion. 



274 D. Discrete wavelet transform and wavelet leaders 

275 Wavelets have been widely used in data analysis 

276 and turbulence research [THHH I23H271 [75H75]. Several 

277 wavelet-based methods have been proposed by several 

278 researchers to extract the scaling exponents from a scal- 

279 ing time series, for example, wavelet coefficients (WC), 

280 WTMM PU [111 [75], WL [231 121 HE], etc. We consider 
28i here WC and WL. 

282 The discrete wavelet transform (DWT) is defined as 



■4>(k,j)= / x(t)<p(2- j t-k)dt 



(24) 



283 where <p is the chosen wavelet, ijj(k,j) is the wavelet co- 

284 efficient, k is the position index, j is the scale index, and 

285 I — V is the corresponding scale [77J [79] . The first way 

286 to detect the scale-invariant properties is to consider the 

287 wavelet coefficients 



z q {j) = MMI 9 > - 2* T <«> 



(25) 



288 where r(q) are the corresponding scaling exponents. 

289 Every discrete wavelet coefficient t^(k,j) can be asso- 

290 ciated with the dyadic interval g(k,j) 



g(k,j) = [tfk,2i{k + l)) 



(26) 



291 Thus the wavelet coefficients can be represented 

292 -0(g) = ip(k,j). Wavelet leaders are defined as 



l(k,j)= sup \ip{g')\ 
e'c3e(k,j),j'<j 



(27) 



293 where 3g(k,j) = g(k-l,j)Ug(k,j)Ug(k+l,j) [2511251175] . 

294 Thus power law behavior is expected 



z,C?) = MM 9 ) - z M9) 



(28) 



295 in which t(q) is the corresponding scaling exponent. Its 

296 efficiency has been shown for various types of data set 

297 [25H2H1EH]. 

298 Let us recall some previous comparison studies be- 

299 tween WTMM, MFDFA and WL. Oswierimka et al. 

300 [69j performed a comparison study between WTMM and 

301 MFDFA by analyzing synthesized data. They stated that 

302 the MFDFA provides a better estimation of singularity 

303 spectrum than WTMM. Jaffard et al. [25] stated that WL 

304 provides a better singularity spectrum than WTMM. Ser- 

305 rano and Figliola [27] performed a comparison study be- 

306 tween MFDFA and WL. They found that WL performs 

307 better than MFDFA. However, for a short time series, 
3oa MFDFA is proposed to extract multifractal spectrum. A 

309 detailed comparison can be found in Ref. [25J, [BS], [27] . 

310 respectively for WTMM and WL, MFDFA and WTMM, 
3ii and WL and MFDFA. 

312 However, we argue here that DWT violates two facts of 

313 the time-frequency representation of a time series. First, 

314 the scale of a time series from complex system, for exam- 

315 pie, turbulent flows, is continuous in a statistical sense, 

316 but not discrete on several scales [3T] [33J . The other one 

317 is that for a certain scale, it may not exist all the time 
3i8 [33] [38j [61] ; see also the discussion in the next section. 

319 Thus to represent a signal by using a DWT is not con- 

320 sistent with the physical aspect. 



III. NONLINEAR EFFECTS 
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FIG. 1. (Color online) (a) a fifth order numerical solution 
(thick solid line) for Duffing equation, (b) An enlarged por- 
tion. For comparison, a sine wave with the same mean fre- 
quency is also shown as a thin solid line. The departure from 
a pure sine wave profile is the result of nonlinear interactions, 
which are nonlinear distortion. 




FIG. 2. (Color online) (a) The first three IMF modes from 
EMD; (b) the corresponding instantaneous frequency u> from 
Hilbert spectral analysis. Note the variation of the instan- 
taneous frequency within one period. It is an intrawave- 
frequency-modulation, which corresponds to a nonlinear in- 
teraction. The instantaneous energy is encoded as a color. 
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FIG. 3. (Color online) Energy spectra (the second-order sta- 
tistical moments) provided by various methods: Fourier anal- 
ysis (x), continuous wavelet transform with db3 wavelet (□), 
first-order DFA (<£>), SF (A) and Hilbert spectral analysis 
(solid line. The wavelet, SF and DFA spectra have been con- 
verted into spectral space by taking / = l/l, f — l/£ and 
/ = 1/n, respectively. Except for the Hilbert-based method, 
all methods require high order harmonic components, which 
are not present in the original signal, see Figs. [I] and [2] to rep- 
resent this nonlinear process. The high-order harmonics may 
lead to an artificial energy flux from low frequencies (large 
scales) to high frequencies (small scales) in spectral space. 



322 We first consider nonlinear effects by using the classical 

323 Duffing equation, which reads 






■x(l 



= bcos(£lt) 



(29) 



324 in which e is a nonlinear parameter. It can be consid- 

325 ered as a pendulum with forcing function 6cos(f2i), in 

326 which its pendulum length varies with the angle. Fig- 

327 ure [l] shows a fifth-order Runge-Kutta numerical solu- 



328 tion (thick solid line) with e = 1, b = 0.1, fl = 2ir/25, 383 gth-order statistical moment) provided by SFs, MFDFA 

329 and [x(0),a/(0)] = [1,1]. The sampling frequency sets as 384 and wavelet are converted from the physical domain into 

330 10 Hz. For comparison, we also show a pure sine wave 385 the spectral domain by taking / = l/£, f — 1/n, and 

331 (thin solid line) with the same mean frequency. One can 386 / = 1//, respectively. The corresponding scaling expo- 

332 see that the wave profile of the solution of the Duffing 387 nents are estimated on the range 0.001 < / < 0.1 (we set 

333 equation is significantly deviating from a sine wave. This 388 here the sampling frequency as 1). Wavelet transform is 

334 deviation is the result of nonlinear interaction, namely 389 performed by using the db3 wavelet. Due to the limit a- 

335 nonlinear distortion [35J I3S]. It is also clearly shown 390 tion of the SF analysis and the HSA, we only consider 

336 that there are no high-order harmonic components in the 391 here the non-negative gth-order moment, q > 0, the left 

337 physical domain [33] . Figure [2] shows the first three IMF 392 part of the singularity spectrum. 

338 modes obtained from EMD decomposition and the cor- 

339 responding instantaneous frequency u from Hilbert spec- 

340 tral analysis. The instantaneous energy is encoded as a 393 A. Fractional Brownian motion 

341 color. The instantaneous frequency u of the first IMF 

342 mode is varying within one period. This corresponds to 394 We have shown in previous works that the arbitrary- 

343 the so-called intrawave-frequency-modulation, which is 395 order HSA can be applied to the fractional Brownian mo- 

344 associated with the nonlinear interactions [3J1 ISHJ [39] . It 396 tion [3U [32]. Here we briefly recall these results. FBm 

345 also clear shows that for a certain frequency, it may not 397 i s a Gaussian self-similar process with a normal distribu- 

346 exist clearly all the time. 398 tion increment, which is characterized by H, the Hurst 

347 Figure |3] shows the normalized energy spectra (or the 399 number < H < 1 [50H53] . Note that the singularity 

348 second-order statistical moments) provided by various 4 oo spectra for the above mentioned methods are 

349 methods: Fourier analysis (x), continuous wavelet trans- 

350 form with db3 wavelet (□), the first-order DFA (<>), 

351 SF (A) and HSA (solid line). For display convenience, a = ('(q), f(a) = mm{aq - ((q) + 1} (30a) 

352 the wavelet, DFA and SF spectra have been converted 

353 from physical domain into frequency domain by taking 401 f or cjFs anc ] 

354 / = 1//, / = 1/n and / = l/£, respectively. We empha- 

355 size here that different wavelet families provide a simi- a = £'(q), f(a) = min{aq — £(q) + 2} (30b) 

356 lar spectral curve (not shown here). As pointed out by q 

357 Huang et al. [55] wavelet transform can be considered as 402 f or ^ ne Hilbert-based method and 

358 an adjustable window Fourier transform. Thus it inherits 

359 the shortcomings of the Fourier transform. We observe a = h 1 (q) - 1, f(a) = min{(a + l)q - h(q) + 1} (30c) 



<-! 



360 that except for the Hilbert spectral analysis, all methods 

36i require high-order harmonic components to represent this <. twa a H 

362 nonlinear process. High-order harmonic components are 

363 not present in the time series (see Figs. [I] and [§. It is a = r '( g ) _ 1; y( a ) = mm {( a + i) g _ T ( q ) + 1} (3Qd) 



364 thus a requirement of the method itself, not the physics 

365 33, 381. This is the main drawback of traditional meth- 



</ 



for WC and WL, respectively. Ideally, we should have 



366 ods, in which the basis is given a priori. Therefore it „ ,,,>_. 

367 is inevitable that one requires high-order harmonic com- TTT r , ,-™ 1. i ri , 1 «ia 1 , 

, ,, ,.„ , ,, , , 406 We performed 500 realizations each oi length 2 lft data 

368 ponents to represent the difference between tnc analyzed . . , . . _, . . TTT ° , 

, , ,. , , , ,i,i., 407 points by applying a lourier- based Wood-Cnan algo- 

369 signal and the given basis. We argue here that high or- ... ~-ji "„°,,. ... . , ,. TT ° , 

. , . , ,, ,.„ . . 408 rithm 84 with H = 1/3, which corresponds to the Hurst 

370 der narmomc components may lead to an artificial energy . ' v, 1 1 , 1 „, , ,, , 

„ r . , ■ n 1 \ i 1 • 1 r 409 number of turbulent velocity. We apply the above men- 

371 flux from low frequencies (large scales) to high Irequen- ,. , ,. , , .. ,. ^ ., . , 

/ ,. , 7 . \ ? m, c ™ tioned methods to each realization of the data series. 

372 cies (small scales m spectral space, therefore, power „. „ , . , ,. ,. . ,. ,. . 
, , , .„.,., . . 11 ,1 ■ 411 1 he final spectra and statistical errors arc then estimated 

373 law behavior, it it exists, may be contaminated by this £ ,. ___ .. ,. „. m , ,. r , % 

. .„ . , n TTT ... , ... . , . 412 from these 500 realizations, ligure 4 show results for (a) 

374 artificial energy flux. We will show this point experimen- „„ ,. ., „ ,« ., _ _ , . _lj — a ( /\\ A 

375 tally by analyzing a temperature data set with strong ' \ , \ y .', ,. . ,, ^ ' i- i- 

376 ramn cliff structures in Sec IVl "" q = 6 ^ V '' ( mlddle ) the corresponding scaling exponents 

376 lailip-Cllll SLILICLUICS ill oeL. V >/\ ,i r\^^n/.i,\n i 

1 — ' 415 ((q) on the range < q < 8, (right) the correspond- 

416 ing singularity spectrum /(a), (b) HSA, (c) DFA, and 

417 (d) wavelet, respectively. The symbols are the same as 

377 IV. VALIDATION AND CALIBRATION m the SF symbo i s . Graphically, all methods provide com- 

419 parable estimation of f{a). However, we note that the 

378 In this section, we will validate the Hilbert-based 420 Hilbert-based method slightly overestimates £(<?) when 

379 method by performing a comparison study of simulated 421 q > 6. Additionally both the first- and second-order DFA 

380 ffim with Hurst number H = 1/3 and synthesized multi- 422 provide slight underestimation of h(q) and seem to prc- 
38i fractal random walk with an intermittent parameter u = 423 diet a systematic underestimation of the Hurst number 
382 0.15. For comparison convenience, spectral curves (or the 424 H . The WC and WL provide almost the same estimation 
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FIG. 4. (Color online) Analysis results of f"Bm with Hurst number H — 1/3. (a) structure functions, left: S q (£) with g = 
(O), q = 2 (□), q — 4 (A) and q — 6 (V); middle: the corresponding scaling exponents £(q) on the range < q < 8; right: 
the corresponding singularity spectrum h(a). (b) results of Hilbert spectral analysis, (c) multifractal detrended fluctuation 
analysis (d) wavelet coefficients and wavelet leaders. The symbols are the same as for structure functions. Scaling exponents 
are estimated in the range —3 < log 10 (/) < —1. The statistical errors are estimated from a total of 500 realizations. 
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FIG. 5. Illustration of the discrete cascade process. Each step 
is associated with a scale ratio of 2. After n steps, the total 
scale ratio is 2™. 



425 for this simple monofractal process. It seems that they 

426 provide a better estimation than Hilbert and DFA meth- 

427 ods. This result is not in full agreement with Oswi^cimka 

428 et al. [SS] , who stated that the MFDFA provides a better 
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FIG. 6. A sample for one realization of length 2 17 points 
with jj, — 0.15, (a) the multifractal measure, and (b) the con- 
structed multifractal nonstationary process. 



429 estimation of H than WTMM. 

430 The above results show that all methods provide com- 432 H = 1/3 . However, it seems that SF and wavelet based 

431 parable prediction of singularity spectra for fBm with 433 methods provide a better estimation. 
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FIG. 7. (Color online) Multifractal random walk with fi = 0.15. (a) Structure function , (b) Hilbert spectral analysis, (c) 
multifractal detrended fluctuation analysis and (d) wavelet coefficients and wavelet leaders. The symbols are the same as in 
Fig. [4] The statistical error bars are estimated from the total f 00 realizations. 



434 B. Multifractal simulation 

435 We show now that the new method applies to multi- 

436 fractal time series. First, let us consider a multiplica- 

437 tive discrete cascade process to simulate a multifractal 

438 measure e(x). Figure [5] illustrates the cascade process 

439 algorithm. The larger scale corresponds to a unique cell 

440 of size L = £0 A", where £q is a fixed scale and Ai > 1 

441 is a dimensional scale ratio. For discrete models, this 

442 ratio is often taken as Ai = 2. The models being dis- 

443 crete, the next scale involved corresponds to Ai cells, 

444 each of size L/\i = ^oA™ -1 . This is iterated and at 

445 step p (1 < p < n) there are X\ cells, each of size 

446 L/\\ — £ \™~ p . There are n cascade steps, and at step 

447 n there are A" cells, each of size £q, which is the smallest 

448 scale of the cascade. To reach this scale, all intermedi- 

449 ate scales have been involved. Finally, at each point the 

450 multifractal measure writes as the product of n cascade 

451 random variables 



stationary multifractal time series as 



e(x) = [J W p , 



(31) 



452 where W PtX is the random variable corresponding to po- 

453 sition x and level p in the cascade [55] • Following mul- 

454 tifractal random walk ideas 86j 87], we generate a non- 



u(x) 



tix'Y^&Bix') 



(32) 
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456 where B(x) is Brownian motion. Taking lognor- 

457 mal statistic for e, the scaling exponent Q{q) such as 

458 {\Au T {t)\ q ) ~ t^ 9 ) can be shown to be written as 



C(«) 



2 { 4 2' 



(33) 



where fi is the intermittency parameter (0 < fi < 1) 
characterizing the lognormal multifractal cascade. 

Synthetic multifractal time series are generated follow- 
ing Eq. (32 1. For each realization, we choose n = 17 lev- 



463 els, corresponding to data sets with data length 131,072 

464 points each. A sample for one realization is shown in 

465 Fig. [6] (a) for the multifractal measure and (b) for the 

466 nonstationary multifractal time series with \x = 0.15. 

467 We perform 100 realizations with intermittent parameter 

468 fi = 0.15. Except for the structure functions, we apply 

469 all methods to each realization by dividing one realiza- 

470 tion into eight subsets with 2 14 data points each. The 

471 spectra for each realization are averaged over these eight 

472 subsets. The final spectra and error bars are respectively 

473 ensemble average and standard deviation estimated from 

474 these 100 realizations. 



475 
476 
477 
478 
479 
480 
481 
482 
483 
484 
485 



489 
490 

491 
492 
493 
494 

495 
496 



Figure [7] shows the results of (a) SF, (b) HSA, (c) 
MFDFA and (d) WC and WL, respectively. The sym- 
bols are the same as in Fig. |4j The theoretical scaling 
exponents and the corresponding singularity spectrum 
f(a) on the range < q < 8 are shown as a solid line in 
the corresponding sub figures. We see that SFs under- 
estimate C(g) when q > 4. The corresponding estimated 
singularity spectrum /(a) deviates from the theoretical 
line when a < 0.4, corresponding to q > 2.5. It also has 
the largest statistical error. Hilbert methodology slightly 
underestimates £(q) when q > 5. It provides a better esti- 
mation of scaling exponents and /(a) than SFs. MFDFA 
provides the smallest statistical errors for spectral curves 
F q {n), scaling exponents h(q) and singularity spectrum 
f(a). However, it still slightly overestimates h(q) when 
q > 6. We note that the first- and second-order DFA 
provide an equivalent result. WC and WL predict al- 
most the same spectral curves, scaling exponents r(q) 
and singularity spectrum /(a). The corresponding sin- 
gularity spectrum significantly deviates from the theo- 
retical curve. We also note that none of these methods 
recover the whole theoretical line on the range < q < 8. 



PASSIVE SCALAR TURBULENCE WITH 
RAMP-CLIFF STRUCTURES 
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FIG. 8. (Color online) A 0.1s portion of temperature data 
showing strong ramp-cliff structures. For comparison, a pure 
sine wave is also shown as a thick solid line. Note that the 
ramp-cliff structure is significantly different from a sine wave, 
which may cause serious artificial energy flux in Fourier spec- 
tral space. 

We now apply the above-mentioned methods to a real 
time data set, a temperature time series as a turbulent 
passive scalar. The data are obtained from a jet experi- 
ment performed at Joseph Fourier University Grenoble, 
France. The bulk Reynolds number is about Re ~ 60000. 
The corresponding Taylor's microscale Reynolds number 
is about Re\ — 250. The initial temperature of the two 
flows are Tj = 27.8°C and T = 14.8°C. The measure- 
ment location is in the mixing layer and close to the noz- 
zle of the jet. The sampling frequency is 50 kHz. The 
total data length is 10 s, corresponding to 500,000 data 
points. Figure [5] shows a 0.1s portion temperature data, 
illustrating strong ramp-cliff structures. For comparison, 
a pure sine wave is also shown. Obviously, the so-called 
ramp-cliff structure is a large-scale structure with a very 
sharp interface [88-91]. We note that the profile of ramp- 
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FIG. 9. (Color online) Energy spectra (or the second-order 
statistical moment) provided by several methods. The inset 
shows the compensated spectra by multiplying the result by 
js/3 f or Hilbert and Fourier, / 8,/3 for DFA and wavelet, and 
f 2 ' 3 for SFs, respectively. For clarity, the curves have been 
vertically shifted. Both Fourier and Hilbert methods predict 
a clear power law on the range 80 < u < 2000 Hz. Due to 
the presence of ramp-cliff structures, the SF analysis fails to 
capture the power law behavior and DFA and wavelet predict 
a short inertial range on the range 100 < / < 1000 Hz. The 
corresponding scaling exponents are j3g = 1.56 for Fourier, 
&(2) = 1.70 for Hilbert, t&(2) = 2.46 for WL with db3 
wavelet, and h$(2) = 2.47 for the first-order DFA, respec- 
tively. 



515 cliff structures is significantly deviating from a sine wave. 

516 Thus for the Fourier-based methodologies, it is inevitable 

517 that one requires high-order harmonic components to 
sis represent their difference, in which the underlying idea 

519 is a linear asymptotic approximation [37] 138] 161] . This 

520 linear asymptotic approximation process thus leads to an 

521 artificial energy flux from low frequencies (large scales) 

522 to higher frequencies (small scales). It means that the 

523 Fourier-based spectrum may be contaminated by this ar- 

524 tificial energy flux. As another direct consequence, the 

525 artificial redistribution of the energy will lead to an un- 

526 real correlation if we consider cross-correlation between 

527 two scales [92] . 

528 The original time series is divided into 122 non- 
529 overlapping segments with 2 12 data points each. The 

530 finally spectra and statistical errors (the standard devi- 

531 ation) are then estimated from these 122 realizations. 

532 Figure [9] shows the energy spectra (or the second or- 

533 der statistical moments) provided by HSA (solid line), 

534 Fourier transform (dashed line), WL (□), the first-order 

535 DFA (o) and SF (0), respectively. The inset shows 

536 the corresponding compensated spectra by multiplying 

537 a Kolmogorov- Obukhov-Corrsin 93 95j nonintcrmittcnt 

538 scaling exponent 5/3 for Hilbert spectrum, Fourier power 

539 spectrum, 8/3 for WL and DFA, and 2/3 for SF, respec- 

540 tively. Except for the SF, all methods display a clear 

541 power law on the range 80 < / < 2000 Hz or 100 < 

542 / < 1000 Hz, a more than one decade inertial range. 

543 The corresponding scaling exponents are £e(2) ~ 1.70 
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FIG. 10. (Color online) Results for passive scalar (temperature) with strong ramp-cliff (a) Hilbert spectral analysis, (b) 
multifractal detrended fluctuation analysis, and (c) wavelet coefficients and wavelet leaders. The symbols are the same as in 
Fig. 
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for Hilbert, /3 e ~ 1.56 for Fourier, r e (2) ~ 2.46 for 
WL and hg(2) ~ 2.47 for DFA, respectively obtained 
by using a least square fitting algorithm. We note that 
only the Hilbert based scaling exponent £e(2) is close 
to the corresponding nonintermittent scaling exponent 
£e(2) = 5/3 [3B]. It is also comparable with the scaling 
exponent of longitudinal velocity in fully developed tur- 
bulence [96] [97] . Due to the presence of strong ramp- 
cliff structures, the SF fails to detect the correct scaling 
behavior. The influence of large energetic structures on 
SF has been studied in detail by Huang et al. [36] and 
Huang [33 . It is interesting to note that DFA and WL 
provide almost the same scaling exponent, which indi- 
cates that the ramp-cliff structure may have the same 
influence on them. We believe that there exists an arti- 
ficial energy flux as we discussed above in both Fourier 
and DFA and WL spectra. Thus they may underestimate 
the scaling exponents [55) . 



562 Figure 10 shows the analysis results of (a) Hilbert spec- 

563 tral analysis, (b) MFDFA, and (c) wavelet transform, 

564 respectively. The symbols are the same as in Fig. |4] 

565 For comparison, the lognormal model of longitudinal ve- 

566 locity [SB] is shown as a solid line in the subfigures of 

567 scaling exponents and singularity spectra. Due to the 



failure of SF analysis, we do not present it here (see 
Rcf. [33 ). Graphically, these three methodologies pre- 
dict power law spectra with small statistical error. The 
corresponding scaling exponents are estimated on the 
range 80 < / < 2000 Hz or 100 </ < 1000 Hz. It 
is found that the corresponding scaling exponents £e(q) 
and singularity spectrum f(a) are close to the lognor- 
mal fitting model, indicating a less intermittent passive 
scalar turbulence field than what was believed before [33] . 
MFDFA and wavelets provide comparable statistical er- 
rors and singularity spectra. Their scaling exponents and 
singularity spectra significantly deviate from lognormal 
model, which is usually considered as evidence that the 
passive scalar turbulence field is much more intermittent 
than the velocity field 88 -"91] . We note that the first- 
and second-order MFDFA provide different scaling expo- 
nents and singularity spectra, which may be associated 
with the different abilities of different order polynomials 

Harm. 



587 As we already mentioned previously, the wavelet and 

588 DFA spectra are strongly influenced by nonlinear large 

589 scale structures (e.g. ramp-cliff structures in passive 

590 scalar turbulence) . Their scaling exponents are thus con- 
59i taminated by high-order harmonics. In other words, the 
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592 statistical property of small scales is contaminated by 647 all the methods require high-order harmonics to represent 

593 nonlinear large-scale structures. We believe here that 648 the ramp-cliff structures. Therefore, the singularity spec- 

594 the scaling exponents and singularity spectrum provided 649 tra provided by DFA and WL are contaminated by this 

595 by them are not correct. Since the HSA has a very lo- eso large nonlinear structure. In fact, it already has been 

596 cal ability in both physical and spectral spaces, together «i reported by several authors that for passive scalar tur- 

597 with the ability of intrawavc-frcquency-modulation for 652 bulence the second-order SF and Fourier power spectrum 

598 nonlinear processes, the effect of ramp-cliff structures is 653 are not consistent with each other [9"Ull9Uiri00| . Warhaft 

599 constrained in the amplitude-frequency space. Therefore, 654 90] stated that "the reason for this is unclear, but appar- 

600 the HSA method may provide a more correct scaling ex- 655 ently stems from the Fourier transform itself." There is 
6oi ponent and singularity spectrum. We note that for DFA 656 no mathematical transform involved in SF analysis. Now, 

602 and the wavelet method, the large deviation from a log- 657 it seems quite clear that not only Fourier-based methods 

603 normal spectrum may be interpreted as a shift problem 658 are strongly influenced by the ramp-cliff structure, but 

604 for moment-based methods when the translational invari- 659 also SF analysis |36j and DFA. 

605 ance is broken [29]. I30j. It seems that the Hilbert-based eeo Our experience is that the HSA is a direct mea- 

606 method can automatically correct this problem. We also mi surement of scale-invariant property in the amplitude- 

607 underline here that the Reynolds number of the present 6 62 frequency space. It thus requires a much larger sample 
60s passive scalar data set is about Rex — 250. Thus the 663 size to get a convergence result than SFs, DFA and WL. 
609 strong ramp-cliff structure may be recognized as an ef- 664 Thus for a sample of small size without large-scale struc- 
6io feet of the finite Reynolds number. We will address this 66 5 tures, SF analysis, DFA or WL are useful to extract seal- 
en issue elsewhere. 666 ing exponents and singularity spectrum since all method- 

667 ologies provide almost the same result. In this situation, 

668 the HSA is useful to catch the scaling trend. However, 
m VI. CONCLUSION 6M if the data set possesses a large-scale structure (e.g., the 

670 ramp-cliff structure in passive scalar turbulence, the sea- 

T , , . ,,. 67i sonal cycle in the daily discharge of river flow, etc.), we 

6i3 In summary, we introduced in this paper a new , . TTrl . . . . ° . . ' 

,, , , ,., , TT .,, , , , ,672 argue that the HSA is the best choice. 

6i4 method, namely arbitrary-order Hubert spectral anal- ° 

6i5 ysis, to characterize scale-invariant properties directly 673 , Fmall y> we , wou } d llk f *« P™ vlde some comments on 

6i6 in the amplitude-frequency space OT31. It is an ex- 674 the moment-based methods, for example, the methods 

6i7 tended version of Hilbert-Huang transform [38l ESI |43] . 675 Presented in this work^nd the singularity-based ap- 

6i8 The main advantage of the Hilbert-based methodology 676 P roadies ^ g ' GMWP J^M)- Due to the convergence 

6i9 is its fully adaptive M\ and very local ability both in 677 P^lem, the moment-based methods seem to require a 

620 spectral and physical domains MM- Thus, it is not 678 much lar S er sam P le size than the singularity-based ap- 

62i necessary to require high-order harmonics to represent 679 Peaches. Thus the statistical error bars presented in 

622 nonlinear and nonstationary processes, which is usually 680 thls P^Per could be non-significant due to the finite sam- 

623 required by conventional Fourier-based methods, such as 681 P le size ' Another issue is the ri S ht P art of the sin S ular - 

624 Fourier transform, wavelet transform, etc. We illustrated 682 ^ s P ectrum > which corresponds to evaluating the neg- 



ative order statistical moments for the moment-based 



625 the nonlinear effect by using the Duffing equation. It is , 

626 found that not only Fourier-based methods, but also SF 6M meth °ds. As we already mentioned previous y, it may 



627 analysis and DFA are influenced by nonlinear processes. 685 £e inaccessible for most of the moment-based methods 

628 It is also found that the HSA can constrain the high-order 686 W- However, the singularity-based methods can over- 

629 harmonics by using the intrawave-frequency-modulation 687 come , this P^Uem [29H30] It seems that the idea of the 

630 mechanism for the nonlinear distortion EH ESI M- 68S ^planty-based approaches can be extended into the 

■t-ir ,, r i • i tii tt-11 689 Hilbert frame. This will be presented in a future work. 

63i We then performed a comparison study of the Hubert- r 

632 based methodology with SF analysis, MFDFA, and WL, 

633 by analyzing fBm simulations with Hurst number H = 

634 1/3 and a synthesized multifractal lognormal random ma ACKNOWLEDGMENTS 

635 walk with intermittent parameter \x = 0.15, respectively. 
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